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Abstract 

Transrectal biopsies under 2D ultrasound (US) control are the current clinical standard for prostate cancer diagnosis. 

The isoechogenic nature of prostate carcinoma makes it necessary to sample the gland systematically, resulting in a 

^^ . low sensitivity. Also, it is difficult for the clinician to follow the sampling protocol accurately under 2D US control 

^N ' and the exact anatomical location of the biopsy cores is unknown after the intervention. Tracking systems for prostate 

biopsies make it possible to generate biopsy distribution maps for intra- and post-interventional quality control and 

3D visualisation of histological results for diagnosis and treatment planning. They can also guide the clinician toward 

non-ultrasound targets. In this paper, a volume-swept 3D US based tracking system for fast and accurate estimation 

of prostate tissue motion is proposed. The entirely image-based system solves the patient motion problem with an 

a priori model of rectal probe kinematics. Prostate deformations are estimated with elastic registration to maximize 

accuracy. The system is robust with only 17 registration failures out of 786 (2%) biopsy volumes acquired from 47 

^^ . patients during biopsy sessions. Accuracy was evaluated to 0.76+0.52 mm using manually segmented fiducials on 

ly^ ' 687 registered volumes stemming from 40 patients. A clinical protocol for assisted biopsy acquisition was designed 

O . and implemented as a biopsy assistance system, which allows to overcome the draw-backs of the standard biopsy 

procedure. 
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C^^ ' 1. Introduction 

o. 

1.1. Prostate biopsies 

Today, prostate biopsies are the only definitive way to confirm a prostate cancer hypothesis. The current clinical 
standard is to perform prostate biopsies under 2D transrectal ultrasound (TRUS) control. The US probe is equipped 
k^ , with a needle guide for transrectal access to the prostate, cf. Fig.[T] The guide aligns the needle trajectory with the end- 

^ ' fire US image plane, which makes it possible to visualize the trajectory on the image for needle placement control. 

However, early- and mid-stage carcinoma are mostly isoechogenic, i.e. not visible in US images, which makes it 
necessary to sample the gland according to a systematic pattern. The standard protocol consists in the acquisition of 
10-12 cores and takes roughly into account that most tumors (about 70%) develop in the peripheral zone of the gland 
(see Fig.|2]l. 

1.2. Clinical issues 

The current standard biopsy procedure has several shortcomings: first, it is difficult for the clinician to reach 
systematic targets accurately because he has to move the probe continuously to place the needle; a constant visual 



reference is lacking 01611 . Second, performing a non-exhaustive systematic search for an invisible target implies 



that the target can be missed. Negative results leave the clinician in a dilemma when the cancer hypothesis cannot be 
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Figure 1: 2D-TRUS guided prostate biopsies. The patient is in dorsal or lateral decubitus and locally anesthetized. The 2D TRUS probe is inserted 
into the rectum to position the needle. A rigidly attached needle guide ensures that the puncture trajectory lies in the US plane. Fig. (b) illustrates 
an 18 Gauge (diameter of 1.27 mm) biopsy needle spring gun. Fig. (c) shows the 2D TRUS image containing the needle (bright line with reflection 
artefacts) and the virtual puncture trajectory known from guide calibration (dotted line). Fig. (a) was found on the web site of the National Cancer 
Institute (www.cancer.gov). Fig. (c) was found on the web site of BK Medical (www.bkmed.com). 
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Figure 2: 12-core protocol. Fig. (a) shows the schematized protocol in the coronal plane of the prostate. Fig. (b) illustrates the protocol in 3D. PZ 
is the peripheral zone, TZ the transition, and CZ the central zone of the prostate. 



discarded: his only option is to repeat the biopsy series. Furthermore, the location of the acquired samples with respect 
to the patient anatomy is only very approximately known after the intervention. Uncertainty about tumor location is 
a major reason why prostate cancer therapy generally consists in the treatment of the entire gland. Incontinence and 
impotence are frequent side-effects of whole gland treatments that considerably reduce the quality of life. For this 
reason a rapidly increasing number of research groups is currently investigating methods that allow the implementation 
of focal therapy strategies [JJ. 

1.3. Proposed solutions 

Several authors have proposed to perform biopsy under MR-guidance to overcome the tumor visibility problem. 



Hata et al il3h . Susil et al. J26] and Krieger et al. JlSf propose an MR compatible end-effector for transperineal biopsy 
core acquisition, while Beyersdorff et al. |6] propose an end-effector for transrectal access. Additionally, Stoianovici 
et al. 1 25] propose a fully actuated transrectal biopsy robot. MR imaging is, however, a costly and sparse resource; it 
is thus unlikely that MR-guidance will become a standard for the millions of prostate biopsies performed in the US 
and the EU alone every year MR-guidance is, however, an interesting option if a biopsy series needs to be repeated 
because of inconclusive results. 



Baumann et al. |3] and Xu et al. 11281] proposed to acquire a US volume before the intervention and to use it as 
an anatomical reference. The stream of US control images acquired during the intervention is then registered with the 
reference volume. This allows the projection of targets defined in the reference volume into the control images, and, 
conversely, the biopsy trajectory, known in control image space, into the reference volume. This technique makes 
it possible to improve biopsy distribution accuracy by showing the current trajectory in a fixed reference together 
with the trajectories of previously acquired biopsy cores. It also allows the clinician to aim for targets defined in the 
reference volume during a planning phase and to obtain the precise biopsy positions for post-interventional analysis. 
An example for non-US targets are suspicious lesions found in MR volumes that are multi-modally registered with 
the US reference volume. It is also possible to derive targets from more sophisticated statistical atlases |24] or, in 
the case of repeated biopsies, they could consist of previously unsampled regions. After the intervention, the biopsy 
trajectories in the reference volume can be combined with the histological results and used for therapy planning. 

Xu et al. [28] acquire a freehand 3D US volume and use 2D control images during the intervention. The control 
images are tracked in operating room space with a magnetic sensor mounted on the probe. In a second step, image- 
based registration is performed in about 15s to compensate for small organ and patient movements. A similar approach 
was proposed by Bax et al. ||5J and Cool et al. I12j] . who use an articulated arm for 2D US beam tracking and 
to acquire a 3D reference volume. Bax does not compensate for patient and gland movements. However, pain- 
related pelvis movements are frequent, since the patient is not under total anesthesia. In this case, both methods 
risk to lose track of the gland because the US beam is followed in operating room space and not in organ space, 
requiring the acquisition of a new reference volume. The acquisition of a 3D volume with a 2D transducer is time- 
consuming and the reconstruction process inevitably introduces some distortion in the volume. Furthermore, the 
US probe continually deforms the gland during needle placement, and it is difficult to estimate these non-linear 
deformations on 2D images. This can lead to inaccurate estimations of the puncture trajectory. Envisioneering Medical 
Technologies commercialize the TargetScan system |2] which uses a side-fire probe, flexible biopsy needles and a 
motorized encoded stepper to place and track the needles. This device shares the draw-backs of the system proposed 
in |l5lll2[]. In yl], we proposed to use a US probe with an articulated, motorized transducer array that allows to acquire 
3D images without moving the probe, which reduces the distortions in the reconstructed volume. The sweep duration 
of 0.5-5 s makes it possible to acquire volumes also during the biopsy core acquisition phase without delaying the 
procedure noticeably. The rich spatial information in these 3D control volumes makes it feasible to design a purely 
image-registration based tracking method that can cope with most types of patient movements during the procedure. 
In ||4|] a deformation estimation step was added to the registration pipeline, which reduced the target registration error 
(TRE) below 1 mm, with a registration time of 6 s. Recently, Karnik et al. il4ll reported a TRE of 6.1 ±2.0 mm 
(RMS) with maximum errors of more than 10 mm for transformations obtained with the articulated arm of the system 
discussed in Ii5ni2] ("pre-registration error"). To achieve a more clinically acceptable error, which they define to be a 
TRE of less than 2.5mm, they evaluated the accuracy of non-linear image-based registration using reconstructed 3D 
TRUS volumes. They obtained a TRE of 1.5±0.8 mm with a computation time of 5-90 s. The proposed registration 
algorithm is based on local optimization and is initialized with transformations obtained from the encoded arm in 



I. Global Rigid 



Transformation 
space 



Model of 

endorectal probe 

kinematics 



DOF 



Image resolution 



Coarsest 



Optimization 
method 



Systematic 
Search 





11. Local Rigid 






\ 


Translation and 
rotation 


\ 






6 








Coarsest 


/ 




/ 


Local/ 
Powell-Brent 





III. Deformation 






\ 


Linear elastic 
deformations 


\ 






-100k 








Coarse to fine 


/ 




/ 


Variational/ 
Full multlgrld 



Figure 3: Registration pipeline. The dimensionality of the transformation space and the image resolution are successively increased. 

operating room space. A new reference volume has to be acquired when the patient moves outside the capture range 
of the registration algorithm. 

In this paper, we present a purely US image registration based prostate tracking system capable of deformation 
estimation. The system does not require beam tracking and reduces hence hardware requirements. It is more accurate 
and faster than other systems that achieve clinically acceptable accuracy as defined in 111 4.1 . and also less sensitive 
to patient motion. The registration framework, previously presented in |l3i|4|], is extended and improved: a method 
is outlined that allows to increase the number of levels in a multi-resolution registration framework by reducing the 
information loss on coarse levels. This makes it possible to reduce registration time significantly. The kinematic model 
||3|], used to estimate plausible positions of the US probe with respect to the gland, is extended to cope with varying 
probe insertion depths. The similarity measure introduced for deformation estimation in |4], capable of dealing 
efficiently with the frequently occurring local intensity shifts caused by probe pressure variations and changing US 
beam angles, is analyzed in more detail. The measure is less complex than the correlation coefficient and therefore 
requires less data to yield statistically robust estimates, which makes it possible to use it on coarser resolution levels. 
In a second part we will present our novel clinical prostate biopsy application which is able to compute precise 
biopsy and cancer maps. It also allows for guidance towards non-US targets and provides the facility to visuahze the 
location of samples acquired during a previous biopsy session. Finally, the different steps of the tracking algorithm 
are validated on a large set of patient data and an in-depth discussion of the clinically acceptable tracking error is 
presented. 



2. 3D ultrasound-based prostate tracking 

2.1. Coarse-to-fine strategies 

Non-linear registration of 3D US image streams is currently the most promising approach to perform organ track- 
ing with deformation estimation. The principal challenges of image-based tracking systems are robust outcomes and 
computational efficiency. A technique to achieve both goals are coarse-to-fine registration strategies that successively 
increase the degrees of freedom (DOF) of the transformation space and the image resolution. An overview of this 
strategy is given in Fig. |3] 

2.2. Coarse-to-fine transformation spaces 

We propose a 3-step registration pipeline (see Fig. [3] upper row) for a purely image-based estimation of the rigid 
and residual elastic transformations between the prostate imaged in a reference volume and in a tracking volume. 
At each step, the number of degrees of freedom of the transformation space is increased. This strategy stems from 
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Figure 4: Geometry of end-fire, volume-swept 3D US view cone. Fig. (a) depicts tlie 2D transducer aiTay plane, (b) tlie motor sweep plane, and (c) 
gives a 3D view of the scanned volume. 

the observation that image registration with high dimensional transformation spaces requires, in general, a fairly good 
initial estimate in order to converge to the physically correct solution. If initial estimates of the required quality are not 
available beforehand, they need to be computed on a lower dimensional transformation space on which the distance 
measure is less exposed to local minima. Ideally, if it is small enough, the search space can be explored systematically, 
in order to find all major local minima. 

In the presented approach, three search spaces are explored using a registration pipehne. The first step of the 
pipeline consists of a systematic exploration of a 3 DOF model of the probe kinematics, that integrates rectal and 
image formation constraints. This model excludes solutions that are not plausible in relation to the prostate. The 
systematic search yields a set of points near strong local minima that are investigated in the second step of the pipeline. 
Here, a local optimization is performed on each minimum using a classic 6 DOF rigid transformation space. Finally, 
the best result is retained as the initial estimate for the third registration step, which estimates the elastic deformation. 

2.3. Coarse-to-fine image resolutions 

The second coarse-to-fine strategy operates on the image resolution. While the first two registration steps are 
executed at very coarse resolution levels of a Gaussian image pyramid, elastic registration descends to finer levels to 
estimate local details like the deformation caused by a needle insertion. Large parts of the transformation are therefore 
computed on very coarse levels, not only boosting the registration speed, but also giving the algorithm a more global 
perspective at early stages of the process. 

Special attention is paid to the construction of the multi-resolution image pyramid. The volume covered by the 
end-fire US beam corresponds to a section of a torus (cf. Fig.HJi. The beam borders in Cartesian space are therefore 
non-trivial and a complex mask is required to differentiate voxels that carry information from background voxels. 
The latter must be ignored during image processing to avoid biases. However, most image processing (e.g. gradient 
computation and Gaussian blurring) is performed on the neighborhood of a voxel, i.e. it cannot be carried out if the 
neighborhood is incomplete. While this is a minor issue for high resolution images, it poses a major problem on 
coarser levels of an image pyramid. The entire border is lost with each down-sampling, and the ratio of border voxels 
to inner voxels increases. As a consequence the information containing volume decreases exponentially with each 
resolution level and the multi-resolution approach would be limited to a small number of levels. 

To cope with this important issue, a first order border extrapolation is performed on the mask borders of a pyramid 
level k before computing level k-\- \. Let / : M-' — > M be the image, p efi^ a. voxel position and ^ e N^* a neighboring 
voxel of p. Finally, let li e N"* be the distance vector d - q- p. Then, the Taylor expansion of I(xd + p),x G R at 
point yields, after solving for I{p), 
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Figure 5: Information loss and image extrapolation. The upper row sliows cuts tlirougli a multi-resolution pyramid of volumes constructed without 
information loss containment, from finest (left) to coarsest (right). Note that the volume of information is significantly reduced on the coarsest level 
(upper row, right-most). The lower row shows a multi-resolution pyramid built with the presented image extrapolation technique. The volume of 
information remains almost constant. 



This expression is discretized by setting x -2, which is the smallest number for which the discrete data prolongation 
is meaningful (it follows that two consecutive voxels in direction d are needed for an estimation). We get, therefore, 



hip) 



2I(p + d)- Up + Id) = 21{q) - I(2q - p). 



(1) 



The final intensity is the mean of Iq(p) for all q e N2b{p) for which the computation is possible, where Noeip) is the 
26-neighborhood of p. If less than a third of the directions can be computed, the voxel is discarded. Fig. |5] shows that 
the volume containing information remains almost constant with this approach. 



3. Rigid registration 

3.1. Parametric optimization framework 

The first two steps of the registration pipeline (cf. Fig. |3]l determine the rigid transformation between the prostate 
in the reference and the tracking images. Both steps operate on low-dimensional search spaces that can be efficiently 
solved with a parametric formulation of the optimization problem. We look for a set of parameters 0* in a parameter 
space c R'' such that 

0'[h,h;ip] = argmin£)[/i,/2 o ip{0)], (2) 



where the /,■ 



are images and ip{6, x) : x 
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defines a spatial transformation in function of Q and 
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. Finally, D[Ii,l2] is an image distance measure. In the presented approach, the Pearson correlation coefficienij 



Drc^ 



Cov(/i,/2) 



VVai-(/i)Var(/2) 
is used as image distance measure for parametric optimization. 
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'Also known as cross correlation or simply correlation coefficient. 
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Figure 6: Model of probe head kinematics under endorectal and prostate image formation constraints. Fig. (a) shows the construction of the origin 
S (0, 0) of the surface discretization, which corresponds to the intersection point of a line defined by the rectal fixed point R and the ellipsoid center 
C with the ellipsoid. Fig. (b) illustrates the exploration of the surface with the probe using the spherical angles a and /J. The US origin O is placed 
on the surface of the ellipsoid. Fig. (c) finally models the rotations around the axis of the US probe. 



3.2. Global rigid registration with endorectal probe kinematics 

The presence of local minima in the lar ge r igid search space make it difficult to use local optimization methods 
(e.g. gradient descent or Powell-Brent, |8l |22|]) directly. The objective of the first registration step is therefore to 
estimate the position of the US probe relative to the organ to find a transformation close to the solution. Instead of 
using magnetic or mechanic tracking systems, a model of the kinematics of the probe under physiological and image 
acquisition constraints was designed. The idea is to profit from the following observations: 

1 . The probe head must be in contact with the thin rectal wall in front of the prostate during image acquisition. 
This is a necessary condition to obtain an image of the gland. 

2. The anal sphincter heavily constrains probe motion in the rectum. 

This makes it possible to define a kinematic model of probe head motion with respect to the prostate surface. 
It is assumed that the center of the probe head lies on the principal axis of the probe (the probe axis). It is further 
assumed that the probe axis always lies on a hypothetical rectal fixed point R eM? which approximates the sphincter 
constraints. Finally, an approximate estimation of the prostate surface in the reference image is needed. The role 
of the surface estimation is comparable to that of a region of interest (ROI), it therefore does not need to be very 
accurate. We thus opted to use a simple ellipsoid, which is manually defined by the clinician on the first volume he 
acquires. This step consists in the placement of an axis-aligned rectangular bounding box around the prostate, which 
can be done with a few mouse clicks. All positions for which the US origin O lies on the surface approximation and 
the probe axis lies on R are considered as plausible probe positions, see also Fig.|6h and b. The contact points can be 
defined with two parameters a and fi using a spherical representation of the ellipsoid. A third parameter A is needed 
to model the probe rotations around the probe axis, see Fig. [6];. 

In the presented approach it is further assumed that the probe head position in the image and the probe head radius 
are known. The precision of the fixed point R is not crucial, and it is possible to use the same fixed point R for all 
patients. It is an approximation that stems from a learning set that was semi-automatically registered: for a given 
patient i of the learning set, /?' was defined as the point with minimal distance to the set of probe axes projected into 
the reference volume after registration. The mean of the R' is an approximation of R that works well in practice (cf. 
Sec.|6] where the accuracy of the model is evaluated). 

The resulting 3D transformation space is small enough for a systematic exploration with an acceptable computa- 
tional burden. To cope with variations in probe pressure, i.e. the insertion depth of the probe into the rectal tissues, 
the model is explored at five dilTerent depths in the probe axis direction at steps of 5 mm. The parameter range is 
configured such that probe tilt angles of up to 30° and arbitrary rotations around the probe axis are considered. The 
systematic search is performed on very coarse resolution levels (4 resolution reductions of factor 2^) to increase the 
performance. 



3.3. Local rigid registration 

The five probe positions with optimal Dec are used as start positions for muhiple rigid registrations. The distance 
measure is optimized on this 6 DOF search space using a fast local method (Powell-Brent |8, 22]). The transformation 
(p{o, a>;x):xt-^ R^x + o models the rigid space, with the origin o e E^ and the three Euler angles w e M^ that define 
the rotation matrix R^^. The origin o is set to the center of the ellipsoidal prostate approximation described in Sec. 13.21 

The goal of this step is to find the local minima of the distance measure Dec in the neighborhood of the initial 
estimates, which makes it possible to make a robust final choice between them. The estimation with minimal distance 
is retained. This step can be executed on a fairly coarse level since the finer details will be considered in the following 
elastic estimation step. The rigid registration can hence be performed rapidly. 

4. Elastic registration 

4.1. Framework for non-linear registration 

Image-based deformation estimation can be formulated as an optimization process of a local distance measure. 
Let /i , /2 : K^ — > M be images, (^ : IR-' — > K.-' the deformation function and the functional D{Ii ,l2',f] a measure of the 
distance between /i and I2 o ip. In contrast to parametric approaches that use basis functions to build the deformation 
function, we will follow a variational approach |27, 19, 10] and define (p(x) - x + u{x), where m : M^ — > R-' is assumed 
to be a diffeomorphism. The deformation could then be estimated by solving the optimization problem 

(fi* = arg min (£[/i , /i; ^]) , (4) 

where the registration energy & simply corresponds to D. Straightforward minimization of a distance measure, how- 
ever, yields poor results in general due to countless local minima, in particular in the presence of noise, partial object 
occlusion and other imperfections in the image data. Unfortunately, US is a particularly noisy modality, which makes 
3D US based deformation estimation vulnerable to local misregistrations. This problem can be addressed by integra- 
tion of a priori models of the expected deformation. This can be done implicitly by adding further energy terms to the 
objective function. In this work, inverse consistency and elastic regularization energies are added. 

4.2. Inverse consistency constraints 

In non-linear image registration, the forward estimation that minimizes S[Ii,l2',tfi] does not generally yield the 
inverse of the backward estimation that minimizes £[/2, /i; i^], i.e. tpoij/ i^ Id with Id :M? —^ M.^, jc i-> x. Introduction 
of Zhang's inverse consistency constraint 1,29.1 

J[^;^]= f Wif, o cp - IdWl, dx (5) 

Jn 

as additional energy penalizes solutions that lead to inconsistent inverse transformations, where Q c M.^ is the regis- 
tration domain in image space. Estimation of the forward and the backward deformations is coupled by an alternating 
iterative optimization 

/+' = argmin(fi[/i,/2;^] + J[^^^]), (6) 

f 

«A*+i = argmin(£[/2,/i;^] + J[/;«A])- (7) 

Concurrent estimation with mutual correction reduces the risk of local misregistrations. 



4.3. Elastic regularization 

The deformation of the prostate caused by probe pressure is mostly elastic, which justifies the introduction of the 
linearized elastic potential [119(1 

r ^ 2 /) 

B[ip] - 6[u + W] = I 7 Y] {dxjUk + dxtU^ + -(div uf dx (8) 

as additional energy, where A and // are the Lame coefficients, which are related to Poisson's coefficient v and Young's 
modulus E via the equations 

(3.1 + 2u)u , A 

^ - and V ^ 



/I + yU 2(/l + ;U) 

Note that the linear elastic potential prevents the estimation of strong deformations when operating with non- 
physical, fractional forces derived from local image dissimilarities. This choice was made to avoid strong misregis- 
trations when I\ and I2 contain local differences for which the image distance metric is not invariant. It is, however, 
possible to use for example the less restrictive fluid regularization to compute complete forces, and to apply linear 
elastic regularization in a second step. The low target registration error (TRE) of rigid registration (cf. Sec.|6]l indicates 
however that the gland rarely gets deformed more than 10%, which is the reason why we did not further investigate 
this approach. 

4.4. Image distance measure 

As for rigid registration, the image distance measure is the driving energy of the deformation estimation process. 
While a global measure is used for rigid registration, deformation estimation requires a local measure to capture the 
transformation details. The degree of locality of a measure is in general determined by the size of its convolution 
kernel. The more local the measure is, the more details can be captured on a given resolution level. In consequence, 
the number of resolution levels that can be used for deformation estimation, and hence the computation time, depends 
heavily on the locality of the measure. On the other hand, the measure must also be able to cope with modality- 
specific image variations like for example intensity changes when the ultrasound beam angle changes with respect 
to an imaged surface. Experiments on patient data have shown that the sum of squared distances (SSD), which has 
a minimal kernel size, is a poor distance measure for deformation estimation in US images. Local intensity changes 
are frequent due to changing US beam angles with respect to the tissues and probe pressure variations. The Pearson 
correlation that was used for rigid registration requires, however, a considerably larger kernel size to achieve statistical 



power, in particular if its complex derivative is used. The generalized correlation ratio introduced by Roche et al. |l23 1 
allows us to derive a more appropriate correlation-based distance measure: 



DcAhJi, /] = p Yi^xix) - f{h{x)))\ (9) 

xeQ. 

where the function / : K — > M models the intensity mapping. When a linear relationship f - al + (3 is, assumed, 
the measure corresponds to the Pearson correlation Dec used for rigid registration (cf. Eqn. 13.1b . Note that Dec 
estimates the optimal a and yS implicitly from the image data. Removing a leads to a correlation measure with only 
one intrinsic DOE, i.e. 

DcR[h,h-J{-)+P] ^T^Yj^h(x)-l2(x)-/3f. (10) 

The parameter/? can be directly estimated on the difference image using a Gaussian convolution ^o-(x). Introducing 
the registration transformation (p and switching to the continuous domain yields 

Dsc[Iul2;ip]:= f (hix) - I2 o ipix) - ih - I2 '^ ip) * gAx)f dx. (11) 

Jn 

This measure requires less data to achieve statistical power compared to the two-dimensional Pearson correlation, 
but it is more general than the identity-assuming SSD with its zero DOE. It is invariant to low-frequency intensity 
shifts between the compared images, which is why we called it shift correlation Dsc- The invariant frequency range 
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Figure 7: Intensity shift correlation model. Fig. (a) shows the floating volume. Fig. (d) the fixed volume (after rigid registration). Fig. (b) the 
elastic registration with the SSD congelation model and Fig. (c) shows the result with the intensity shift congelation model and inverse consistency 
constraints. All other parameters were identical for both registrations. The SSD driven registration is incoiTect because of various local intensity 
shifts that are caused by the difference in probe pressure between the acquisition of the fixed and the floating images. The intensity shift conelation 
model correctly handles this problem and converges to the physically correct result. 



is controlled by the standard deviation cr of Q. If (j gets smaller, the cropped range gets larger, and the registration 
convergence rate decreases and may even stall if only high frequency noise is left. When used with a multi -resolution 
solver on a Gaussian pyramid (cf. Sec. 14. 51 ). which implicitly performs a low -pass filtering of the intensity variations on 
coarse resolutions, this approach transforms to a band-pass filtering on varying frequency bands. In this configuration 
it is sufficient to chose relatively small standard deviations, without risking registration inefficiencies. Fig.Qillustrates 
the performance of the shift correlation model combined with the inverse consistency constraint. 

4.5. Solver 

Combination of the energy terms yields the alternating system 



f 
ij/* = dxgmm{Dsc[hJ\,>lj]+&[>lj]+IW;<lj])- 



(12) 
(13) 



An iterative two-step minimization scheme is used to solve both objective functions. The Euler-Lagrange equations 
of Eqn.[T2]and[T3]are rewritten as a fixed point iteration 



Jt 

Jt 



i:[/] +. fee [A, /2; /]+//[/; ^'], 
^ X[/] + /©, J/2, /i ; /] + /j[/; ^'l, 



where f e M controls the discretization granularity, and with the elliptic partial differential operator 

■C[ip] - X.[u + Id] - iiAu + {A+ ;[/)Vdiv u. 



(14) 
(15) 

(16) 



which is obtained from the Gateaux-derivative of & at (f, cf. 1 19]. The Gateaux derivative of J at ^ yields the force 
term 

/j[^^;^]=(V(A)o^-(iA°^-W). (17) 

Finally, the Gateaux derivative of Dsc at (f provides us the image-based forces 



h o ip + 2I2 o (fi * go- - h ° IP * Go- * Go-) ■ i^h) ° f- 



(18) 



10 



To avoid the expensive double convolution we use the approximation 

h*0a*0a^h*0a, (19) 

which simplifies the gradient to 

fDscUuh^f] ^-{h-h*Qa-h°V + hoV*&o-)- (V/2) o if. (20) 

Note that ( |20] | can also be obtained by interpreting p as an independent function fi : M? ^> R, and by applying the 
SSD distance measure on the image pair {I\{x), l2i(p{x)) +/3{x)}. The gradient obtained from the Gateaux derivative of 
this measure at {(p,/3), 

f = -{h-l2^V+/3)-l^^^-^^°A (21) 

can be used for an alternating estimation of tp andyS. If /? is estimated from the image instead of computing it with 
the gradient, i.e. by setting P(x) - I2* QAx) — I\ o ip * Qa-{x), we get again ( l20l l. The approximation ( fT9] l can hence 
be analytically justified by using an alternative formulation of the similarity measure (see jlOll . p. 98, for a similar 
discussion on the CC). 

Finally, an iterative algorithm is used to estimate the displacement fields: 

1: while not converged do 

2: compute foscVuh; /] and //[i^*; /] 

3: compute fvsclh, tuip''] and /j[<^*; lA*] 

4: solve Eqn. [T4]for 1^*^' 

5: solve Eqn. [T5]for (/r^^' 

6: end while 
The forces are considered as constants for the resolution of the PDEs [T4l and [TSl and the forward and the backward 
estimation mutually correct themselves at each force update. The PDEs are solved using Red-Black Gauss-Seidel re- 
laxation II9II22I] and the full multigrid strategy is used to obtain a linear computational complexity [7,^. Convergence 



is achieved if the difference of the La-norm of the sum of all local forces is below a threshold for both the forward and 
the backward estimation between two iterations. Note that the algorithms iterate until convergence at each level, de- 
viating thus from the classical multigrid scheme. This is necessary since relaxation is performed on fractional forces. 
Fixed edges and bending side walls are used as boundary conditions on the registration domain lll9|] . 

4.6. Parameter settings 

A recurring problem in variational image registration frameworks is the proper scaling of the energy terms. In the 
presented approach, line searches in the inverse gradient direction are performed within a limited perimeter in order to 
obtain locally exact solutions for the image distance and the inverse consistency forces. The forces obtained are then 
summed together, i.e. they are equally weighted but not averaged. The compressibility constraint of the linear elastic 
tensor poses some problems in the presented framework: since we are operating with fractional forces, it behaves 
like a strong shape constraint and registration will stall when it is enabled from the beginning. We therefore allow 
the deformation estimation to converge with the compressibility constraint disabled (Poisson's v - 0), and turn it on 
only in a second run where local over-stretching and squeezing is corrected (v - 0.48). Young's modulus E has no 
physical meaning in this approach since fractional displacements are used as artificial forces. We set both E and the 
PDF discretization 6t to 0.5. The final displacement <5i^(x)*^' that is added to (p{x)'' is finally limited to ||5i^(x)|| < 1 
to ensure that no jumps over intensity barriers are possible. This constraint greatly improves the overall numerical 
stability. 

5. CUnical appUcation 

In the previous sections we proposed a tracking method for prostate motion using 3D TRUS. In this section a 
clinical application for prostate biopsy tracking and guidance, based on this method, will be presented. The primary 
objective of this application is to combine clinical accuracy with an intuitive user interface and a lightweight clinical 
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Figure 8: Application prototype. The figure shows a prototype of the biopsy assistance application during a biopsy session. 




Figure 9: Clinical flow chart of the biopsy application. Biopsy planning and acquisition phase. 



work-flow. In particular, care was taken such that no logistical and interventional overhead is added to the current 
standard clinical procedure. The application aims to provide solutions for the clinical issues of TRUS prostate biopsies 
enumerated in ll.2l This includes the ability to target suspicious lesions identified on MR images, to visualize locations 
sampled during a previous biopsy session, and to provide interventional biopsy maps as well as post-interventional 
cancer maps. Fig. |9] gives a summary of the clinical work-flow for biopsy acquisition, while Fig. [TOl describes the 
post-biopsy work-flow from histological analysis to therapy planning. The different steps are described in detail in 
the following sections. An illustration of the clinical setup is given in Fig. [8] 



5.1. Panoramic volume as anatomical reference 

The first part of the clinical protocol is performed a few minutes prior to biopsy collection, and consists in the 
acquisition of the anatomical reference. Unfortunately, a single 3D volume of the prostate does not typically contain 
the entire gland because of the pyramidal form of the 3D US beam that cuts the gland at its lateral borders and, in 
particular, near the probe head where most tissue samples are taken. For this reason, three volumes are acquired from 
different positions, elastically registered and compounded into what we name a panorama volume. The bounding 
ellipsoid required for the kinematic model is defined by the clinician after the acquisition of the first volume. 
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Figure 10: Clinical flow chart of the biopsy application. Histological analysis, diagnosis and therapy planning phase. 




Figure 11: Planning phase. Fig. (a) shows the transversal and (b) the coronal semi-automatic segmentation result. Fig. (c) shows the 3D repre- 
sentation of the segmentation including the manually segmented points (yellow). Fig. (d) finally shows the result of non-linear surface registration 
between the MR image and the US reference. 



5.2. Biopsy planning 

The second part of the protocol consists in the definition of non-ultrasound biopsy targets in the anatomical refer- 
ence. This step needs to be simple and rapid since it is not convenient for the clinician to interact with the computer 
while holding a probe, and the patient discomfort increases with the duration of the rectal penetration. The clinician 
therefore first identifies the targets, then he inserts the probe and acquires the reference volume, and finally the soft- 
ware registers the targets with the reference volume. If the target registration requires user interaction, the probe can 
be held in place with an articulated arm in order to free the physician's hands. 

In the concrete case of MR lesion targeting, the segmentation of the prostate in the US reference, required for 
MR-US registration, needs to be performed rapidly for the same reasons. We use the algorithm proposed by Martin 
et al. in [17] to perform MR to 3D TRUS prostate volume fusion. This algorithm uses the MR image prostate 
segmentation as shape prior for a fully automatic detection of the prostatic capsule in 3D US volumes. The MR shape 
is automatically obtained with a probabilistic atlas and a spatially constrained deformable model, using the method 
presented by Martin et al. in JIS ll. cf. Fig.[TT] After MR-US fusion, the segmented MR target can be projected into 
the anatomical reference of the tracking system. The application also provides the facility to project biopsy maps into 
an MR image, using the same techniques. 

In the case of the repetition of a biopsy series, it is interesting to know the exact location of previously acquired 
samples. This is performed by registering the anatomical reference of the previous intervention with the newly ac- 
quired reference using the elastic tracking algorithm presented in Sec.|2] However, we do not yet have enough data 
on repeated biopsies to prove the viability of this approach. The time lapse between two series can be important, and 
the organ can be altered in between, for instance because of tumor growth. Also, the involved imaging hardware can 
be different or differently configured, which could yield dissimilar images. If purely image-based registration should 
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(a) (b) (c) 

Figure 12: Biopsy maps. Fig. (a)-(c) sliow different views of a biopsy map with color-coded Gleason score. 



lack stability, a mixed image and surface-based registration will most certainly be sufficient. 

Finally, it could be interesting to indicate targets with a high probability of carcinoma presence. This can be 
performed by registering a statistical cancer atlas, such as that developed by Shen et al. 1.24.] . with the prostate shape 
segmented in the reference volume. This could lead to a better statistical sampling and could potentially reduce the 
number of biopsy acquisitions per session, see also the optimal protocol for transperineal biopsies developed by Shen. 

5.3. Biopsy acquisition 

After registration and projection of the targets, the biopsy core acquisition phase is started. Current 3D US 
scanners can acquire 0.5-5 volumes per second, depending on the resolution. However, continuous 3D volume streams 
cannot be processed at the same frame-rate with the current implementation of the prostate tracking algorithm. Also, 
the higher frame-rates of 2D US are visually more pleasing. As a consequence, 2D US is used for needle placement. 
Single 3D volumes are acquired only when positional information about the targeted or collected sample is required. 
Volume acquisition can be triggered out of 2D US mode with a foot pedal. The clinician's hands are free for tasks 
like holding the probe and the biopsy gun. The volumes are automatically transferred to the application via TCP/IP 
and DICOM and after reception they are automatically registered. The puncture path that corresponds to the probe 
position when the volume was acquired is then projected into the reference volume and visualized with the planned 
targets. This makes it possible to validate and correct the puncture trajectory before firing the biopsy gun. 

When the needle is correctly placed, the tissue sample is collected by triggering the spring needle gun. Before 
removing the needle, the clinician acquires a US volume containing the needle image. The needle is then automatically 
segmented, if present. It is projected into the reference anatomy, giving the clinician an immediate feed-back about 
the sampling position with respect to previous biopsies. The biopsy distribution can hence be assessed during the 
intervention and additional samples can be collected, if necessary. A typical interventional biopsy map is given in 
Fig. [HI 

5.4. Multi-modal biopsy and cancer maps 

After acquisition, the samples are sent to the pathologist for histological examination. The Gleason score, a visual 
tissue grading scale that is correlated with the aggressiveness of the carcinoma, is determined for each sample and 
entered into the biopsy application. A color code is used to visualize the cancer grade of a sample in the biopsy maps. 
This provides an instant overview over the cancer distribution for diagnosis and treatment planning. 

With the presented system it is possible to localize the biopsy samples accurately. This capacity could help 
to improve existing treatment strategies and might be useful for experimental focal therapy approaches. It can be 
interesting, for example, to project the cancer map onto a MR planning volume. This can be done using the same 
approach as the MR target projection onto the US reference. 
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Figure 13: Viewer for visual validation. The viewer shows transverse (a), sagittal (b) and coronal cuts (c) simultaneously. The registration volumes 
are overlaid: the upper-left and lower-right comers of each view show the tracking volume, while the lower-left and upper-right comers show the 
reference volume. The green cut point is the 3D point where the three image planes intersect. The user can move the cut point freely in each view, 
allowing him to explore the entire volume (changing the point in one view changes the spatial position of the other two views). 

6. Experiments and results 

The tracking system was tested on real patient images that were acquired during biopsy sessions at the Pitie 
Salpetriere hospital, Paris, France. The patient study was approved by the ethical committee of the Pitie Salpetriere 
and was performed with the consent of each patient. The volumes were acquired with a GE Voluson equipped with an 
endorectal RIC5-9 probe. Three volumes were acquired a couple of minutes before biopsy acquisition for the creation 
of a panorama volume that contains the entire gland (see lS.ll l. Interventional volumes were acquired after each biopsy 
shot with the needle left inside. Besides these points, the standard clinical protocol was not altered. Registrations were 
performed on a standard PC with 4GB of RAM, 3GHz processor frequency and two cores. The acquired volumes had 
a resolution of 199^ voxels. A 5-level multiresolution pyramid was used for registration, the resolution on the coarsest 
level being 13^ voxels. 



6.1. Kinematic model 

The first experiment evaluates the capacity of the kinematic probe movement model to estimate a probe position 
sufficiently close to its real position with respect to the gland. For this test, 786 registrations of the biopsy volumes of 
47 patients were performed with the corresponding panorama volumes. The registrations were then visually validated 
either by the clinician immediately after biopsy acquisition or off'-line by one of the authors. This was carried out 
using a volume viewer that allows to overlay and to explore the reference and the tracking volume after application 
of the registration transformation (see Fig. [T3] l. Note that the accuracy of a large and randomly chosen subset of the 
registrations classified as valid was evaluated to 0.8+0.5 mm, cf. the accuracy study presented in Sec. 16.21 which 
indicates that only few registrations were falsely identified as correct. In this study, 769 (97.8%) volumes were 
classified as valid, and 17 (2.2%) were classified as failures. The 17 failures occurred with volumes that did not 
contain enough information about the prostate, i.e. the tracked object was literally "out of view". This was caused 
by inadequate US depth or probe pressure (prostate capsule not visible) in 11 cases, partial probe contact with the 
rectum in 1 case, extremely lateral volumes containing only a small part of the prostate in 4 cases and an incomplete 
panorama in 1 case. Note that the failures were not caused by patient movements. 

A second experiment studies the role and performance of the kinematic model (pipeline step 1 in Fig. [3]) with 
respect to local rigid registration (pipeline step 2). We therefore compute the angular and the Euclidean distances of 
the best transformation predicted by the model from the transformation estimated with local rigid registration to see 
how accurate predictions of step 1 are and how much the second step improves these predictions. The best prediction 
of the model is defined as the transformation that is closest to the transformation produced by step 2. The distance 
between two transformations is defined as the root mean square (RMS) Euclidean distance of the 6 intersection points 
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(a) 



(b) 



(c) 



Figure 14: Probe positions. The Fig. (a) to (c) show the probe positions during biopsy acquisition for three different patients, projected into 
panorama space after registration. The interface between blue and yellow tubes corresponds to the position of the transducer origin, the yellow tips 
correspond to the contact point of the probe with the rectum. The green sphere is the a priori fixed point. One can see that the angular error of the 
probe position predicted by the kinematic model stems from the fact that the fixed point model of the rectal sphincter is rather approximate. 

of the probe head with the canonical coordinate axis, centered at the transducer origin, after respective application 
of the transformations. The study was performed on 269 successful volume registrations stemming from 14 biopsy 
sessions. The ellipsoid surface was discretized using a 16x16 grid, and the 360° rotational space around the probe 
axis was discretized using 24 steps. The model was evaluated at five different depths. The mean Euclidean distance 
of the transducer array origin was 2.0+1.5 mm, and the mean angular distance was 9.7+5.1°. Concerning the five 
different depths at which the model was evaluated, 15% of the best predictions were found at an offset of -10 mm, 
23% at -5 mm, 35% at mm, 21% at +5 mm and 6% at +10 mm. The model hence predicts the position of the 
transducer array quite precisely, which indicates that it is relatively robust with respect to the accuracy of the user- 
defined ellipsoid. Note that evaluating multiple depths can compensate ellipsoid placement errors in the probe axis 
direction. The angular predictions are less accurate, which stems from the resolution of the rotational space around 
the probe axis of 15° and from the fact that using an a priori fixed point as a model for the constraints exerted by the 
rectal sphincter is rather simplistic. After rigid registration, the average distance of the probe axis from the a priori 
fixed point was 4.3+3.1 mm. Fig. [T4l shows typical probe positions during biopsies with respect to the fixed point 
for three diff'erent patients. In conclusion, the results given in this paragraph reflect the trade-off" between accuracy 
and speed of the kinematic model. The role of the kinematic model is to initialize local rigid registration efficiently. 
However, to achieve clinically acceptable accuracy, local registration is mandatory. 

6.2. Target registration error 

This third study evaluates the accuracy of 687 registrations stemming from 40 patients that were classified as 
correct. Fiducials were first manually segmented in the panorama volumes, and then in the registered tracking volumes 
if they were visible. Finding unambiguously identifiable fiducials in the panorama and the majority of the tracking 
volumes was a challenging task. Segmentation consisted in the definition of the barycenter of a structure. Small 
structures like calcifications or kysts were preferred since their limited volume facilitates manual barycenter definition 
and hence reduces the fiducial localization error (FLE), see Fig. [15] We did not segment (parts of) the prostatic 
capsule since it is difficult to identify anatomically corresponding points on surfaces in US images, i.e. the FLE of the 
point distance measure would have increased. Surface-distance measures on the other hand underestimate the tissue- 
correspondence error because they are insensitive to on-surface misalignments. In total, 147 reference fiducials were 
segmented in the 40 panorama volumes (3.7 fiducials per volume), and 1889 corresponding fiducials were segmented 
in the 687 tracking volumes (2.7 fiducials per volume). The statistical power of the simple descriptive statistics that 
were used in this study is ensured by the total number of 1889 evaluated samples, i.e. the small number of samples per 
registration pair is compensated by the large number of evaluated registrations. The FLE was estimated to 0.35+0. 19 
mm via multiple segmentations of the reference fiducials. The anatomical distribution of the reference fiducials is 
illustrated in Fig. [16] It is important to note that the fiducial segmentations are not used by the registration algorithm; 
the fiducial registration error is hence a valid estimator of the target registration error (TRE). 
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Figure 15: Fiducials. Fig. (a) and (b) show corresponding calcifications in the tracking and the panorama volume, Fig. (c) and (d) illustrate a kyst. 







mean 


RMS 


mean dist. 


execution 






distance 


distance 


(worst decile) 


time 






[mm] 


[mm] 


[mm] 


[s] 


1 


unregistered 


13.8±7.9 


17.0 


17.9 


- 


2 


rigid 


1.4±0.8 


1.5 


3.0 


2.1 


3 


deformation 


0.8±0.5 


0.9 


2.0 


6.8 


4 


def. w/o inv. cons. 


1.0±0.7 


1.2 


2.7 


6.4 


5 


def. w std SSD 


1.6±1.5 


2.2 


5.2 


4.3 



Table 1: Accuracy study. Line 2) and 3) present TRE and registration time of the rigid and deformation estimations. Line 4) shows the result for 
the deformation estimation without inverse consistency constraints, and line 5) for the deformation estimation using standard SSD instead of shift 
correlation as similarity measure. The mean en'or in the 'worst decile' was computed on the 10% of fiducials with the largest TRE. 



To evaluate the TRE, the distance between the centers of corresponding fiducials were measured using the Eu- 
clidean distance 



e{P,k)^\\Fl{k)-^f{hJ„\oFl 



(22) 



where P is the patient, F^{k) is the A;th fiducial in the panorama image /q and Ff (k) is the anatomically corresponding 
fiducial in the tracking image /„,. 

The results for both rigid and elastic registration are given in Tab.[T] Rigid registration yields a TRE of 1.4±0.8 
mm while the TRE after deformation estimation is 0.8+0.5 mm. The mean error in the worst decile is 3.0 mm 
for rigid registration, which is reduced to 2.0 mm after deformation estimation. The average computation time of 
rigid registration is about 2 s versus about 7 s for deformation estimation. We also evaluated the performance of 
the deformation estimation without inverse consistency constraints (TRE of 1 .0+0.7 mm) and when using standard 
SSD instead of the proposed intensity-shift compensating metric Dsc (TRE of 1.6±1.5 mm), everything else being 
identical. Both techniques improve thus accuracy. SSD-based registration yields results worse than rigid registration, 
which proves that this measure is inadequate for deformation estimation in US images. On the other hand, shift 
correlation significantly improves the rigid registration result while its computation time is still reasonable. It is hence 
an interesting alternative to more complex similarity measures. 

6.3. Rigid vs. elastic registration 

A comparison of the mean decile TRE of rigid and elastic registration was carried out to evaluate the clinical 
relevance of deformation estimation, see Fig. [T7| Elastic registration improves TRE in all deciles about at least one 
third. In absolute values, it yields the strongest benefit in the worst decile, where the average error is decreased by 1 
mm to 2 mm. The error curve of rigid registration indicates that the gland is barely deformed in the majority of the 
volumes that were analyzed. However, in about 20% to 30% of the volumes stronger deformations can be observed 
that can be reduced with elastic registration. The cross correlation between the TRE and the fiducial distance from 
the transducer is -0.09 for rigid registration, i.e. these two variables are uncorrelated. The cross correlation between 
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Figure 16: Fiducial distribution. This figure shows the distribution of the 147 fiducials that were segmented in the panorama volumes in normalized 
space (unit circle). Normalization was performed by scaling the bounding ellipsoids on the unit sphere, followed by a projection of the fiducials 
onto the transverse mid-plane of the sphere. 

the TRE and the fiducial distance from the prostate center is 0.17 for rigid registration, i.e. the error at the capsule is 
of similar extent than the error at the center. The cross correlation between the TRE and the fiducial distance to the 
needle tip is 0.02, and between the TRE and the fiducial distance to the entire needle trajectory is 0.0. We therefore did 
not measure significantly stronger errors near the needle than elsewhere. This observation has to be, however, taken 
with care since it is possible that including all fiducials to measure a local phenomenon (the volume of the needle is 
very small compared to the volume of the prostate), together with the sparse number of fiducials per volume, might 
compromise the statistical power of this particular test. A visual illustration of the deformation estimation is given in 
Fig.im 

7. Discussion 

We have presented a 3D TRUS based tracking system for prostate motion that occurs during transrectal biopsy 
acquisition. The presented approach differs from existing systems in that it uses a motorized US transducer for 3D 
volume acquisition and in that it is purely image-registration based. A kinematic model of the rectal probe motion is 
used to find the position of the probe with respect to the gland, and not in operating room coordinates, as it is the case 
with approaches that use probe tracking hardware. A registration framework capable of estimating the rigid and the 
residual elastic motion of the prostate was presented. A simple but effective local similarity measure was introduced 
to drive the elastic registration process as efficiently as possible. Special extrapolation techniques were proposed to 
cope with information loss problems in multi-resolution image registration. A clinical application was built on the 
top of the tracking system, allowing the computation of precise biopsy and cancer maps, MR-TRUS target projection, 
MR biopsy maps for therapy planning, and rudimentary support for guidance towards non-ultrasound targets. In this 
section, we will discuss the proposed system. 

7.1. Target registration error 

A prostate tumor is considered as being clinically significant if it has a volume of at least 0.5 cc |l]J,|2lJ]. If the 
shape of the tumor resembles a sphere, this corresponds to a radius of roughly 5 mm. Let us further assume that 
the error is normally distributed. Based on these assumptions, Karnik et al. claim that a tracking TRE of 2.5 mm 
(RMS) yields a probability of 95.4% that the registered targets will lie inside the clinically significant 5 mm radius 



Ill4ll . However, this is only valid if there are no other sources of error in the clinical application. 
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Figure 17: Rigid vs. Elastic. The Fig. shows the error deciles after rigid and elastic registration. 




(a) 



(b) 



(c) 



(d) 



Figure 18: Rigid vs. Elastic. Fig. (a) and (c) show typical rigid registration results, the left half of each image showing the tracking image, and the 
right half showing the reference. In both cases, different probe pressure was applied when the volumes were acquired. The non-linear compression 
cannot be corrected by the rigid registration and local mismatches subsist in particular near the probe head. Fig. (b) and (d) show the result of the 
deformation estimation, which connects the local errors. 
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If the clinician wants to be guided towards non-US targets, for example suspicious lesions in MR images, he has 
to segment and register them with the US reference volume. We can assume that MR target segmentation, MR-TRUS 
registration and TRUS-TRUS registration are statistically independent. In that case, the total RMS error is 



-JL' 



ef. (23) 

If we assume that manual MR segmentation is performed with an RMS error of 1 mm, MR-TRUS registration with an 
error of 2.5 mm 1 17, 18] and TRUS-TRUS registration also with 2.5 mm, we get a total RMS error of 3.7 mm. This 
corresponds to a 82% probability to hit the target. 

However, we also have needle deflection, guide calibration errors, needle depth tracking errors, tissue deformation 
during needle insertion and the error of manual reproduction of the proposed trajectory. Unfortunately we do not 
know the extent of these errors. Let us optimistically assume that each of these factors add an additional RMS error 
of 1 mm. Then we get a total RMS error of 4.3 mm or a probability of about 75% to hit the target. This means that the 
sensitivity of MR-TRUS guided biopsies would correspond approximately to the sensitivity of a systematic sampling 
of the gland. With this error we would thus not gain anything. 

It is therefore necessary to minimize the TRE of each part of the processing chain wherever possible. Reusing the 
assumptions of the previous paragraphs, the total RMS error of the presented system is 3.6 mm (RMS, 83% probability 
to hit the target) with elastic registration, and 3.8 mm (RMS, 81% probability to hit the target) with rigid registration. 
With elastic registration, the probability of missing the target is reduced by 25% compared to a tracking system with a 
TRE of 2.5 mm (RMS). If the TRE of the presented system could be reduced to 0.25 mm (RMS), which corresponds 
to the error of optical or magnetic tracking systems, the probability to hit the target would increase to 85%. In other 
words, further reducing the tracking error by a factor of three improves the hit rate by only 2.5%. For MR targeting, 
future efforts should thus be concentrated on the reduction of the MR-TRUS fusion error. 

Note that the error analysis is not the same for biopsy maps, where the total error is composed only of the TRUS- 
TRUS registration error and the needle segmentation error If we assume, again, that a manual segmentation can be 
performed with an RMS error of 1 mm, the total error would be 1 .34 mm (RMS). If the biopsy map is used for therapy 
planning, however, and if it is projected onto the therapy planning volume, the projection error adds to the total error. 

In conclusion, to achieve a clinically satisfying total TRE, all sources of error must be minimized further. With 
the presented system, the most important sources of error seem to be MR-TRUS fusion and the human error, i.e. 
the capacity of the operator to reproduce a suggested trajectory. Future effbrts should hence be concentrated on the 
optimization of these errors. 

7.2. Rigid vs. elastic tracking 

Elastic tracking is the preferred solution concerning the minimization of the total TRE of the system. However, it 
requires a three to four fold increase in computation time to reduce the risk of missing an MR target, estimated with 
the assumptions in Sec. 17.11 by only about ten percent. Other sources of error are predominant in the system. For 
biopsy maps, where the computation time is not an issue, and for the registration of the panoramas of two different 
sessions, elastic registration is the preferred solution. 

It is interesting to have a look at the error distribution of rigid and elastic registration (cf. Fig. \T7\ . Elastic 
registration reduces the rigid error by almost 50% in the deciles 1-9. In the worst decile, the error reduction is 30%. 
The error is reduced significantly in all deciles, which shows that the proposed non-linear registration method is 
robust. The rigid registration result is not degraded. 

7.3. Patient motion 

A considerable advantage of the proposed approach is that it is relatively robust with respect to patient movements. 
Some patients feel pain during biopsy acquisition, in which case hip displacements due to muscular contractions can 
occur. Other patients are not comfortably positioned and change their pose during the procedure. In the presented 
experiments, registration failures were not due to patient movements, but to a lack of prostate information in the 
tracking volume, which is comparable to a marker that is out of view of an optical tracking camera or to a magnetic 
sensor that is outside the magnetic field. This contrasts with systems that use probe tracking to initialize image-based 
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registration, for which the patient/prostate can move outside the capture range of the registration algorithm. If this 
happens, these systems have to acquire a new anatomical reference volume, which is time consuming and requires 
transitive registration of the previous results. A nice side-effect of a purely software-based and free-hand tracking 
system is that less encumbering hardware is required in the operating room. 

Patient movements can, however, still pose problems with the presented system when they occur during 3D TRUS 
acquisition, which leads to distorted volumes. The same applies when the operator moves the probe during volume 
acquisition. Patient and probe movements can also occur during registration, which is problematic when the system 
guides the clinician toward a target. This can be improved by reducing the registration time, cf. also the discussion in 
Sec.lT4l 



7.4. 3D TRUS vs. 2D TRUS based tracking 

3D TRUS and 2D TRUS based prostate biopsy systems have both strengths and weaknesses. A large incon- 
venience of 3D US and deformation estimation is the volume acquisition time of 0.5-5 s (depending on the image 
quality) and the registration time of currently 7-8 s. While this is still sufficient to give the clinician an intra-operative 
feed-back, it is difficult to implement a system that can guide the clinicians to targets. It is conceptually straight- 
forward to parallelize the registration algorithm on specialized hardware like modern scalar processing units. In this 
case, an increase in speed by a factor of 20-30 seems to be realistic. However, real-time TRUS volume acquisition is 
currently a blocking issue for conveniently fast updates of the tracking position. 

2D TRUS based tracking is less accurate and more sensitive to patient movements. When using 2D TRUS with 
image-based registration, there is a risk that lateral biopsies will not be correctly registered since only a small part 
of the image will contain the gland. The capture range of such a system is smaller than the range of a 3D TRUS 
based system. Another inconvenience of 2D TRUS is the less accurate and more time-consuming acquisition of the 
initial 3D volume that serves as anatomical reference. The impact of the volume reconstruction error cannot be fully 
evaluated with point fiducial based TRE evaluations. Phantom studies have shown distance errors of 3%-4% 151 1 1 III : 



on real patients the accuracy is unknown. Perhaps a combination of 2D and 3D TRUS tracking would yield the best 
of both worlds. 

8. Conclusion 

The biopsy application that we implemented for clinical use gives the clinician an immediate and precise feed- 
back about the current biopsy distribution. The biopsy maps can be fused with the results of the histological analysis 
for diagnostics and therapy planning. The maps of a previous biopsy session can be projected on the current panorama 
image. Finally, the system enables guidance towards targets segmented on MR images. The system only minimally 
alters the current standard procedure and does not add significant logistical complexity to the intervention. 

The cancer map module provides significant clinical value for diagnosis and therapy planning. It makes it possible 
to implement experimental focal therapy strategies based on the localization of positive samples. It is for example 
feasible to perform additional targeted biopsies in a second session around a positive sample to get a very precise idea 
of the cancer distribution. Inversely it is possible to confirm the absence of cancer locally with supplementary targeted 
biopsies. The possibility to project the cancer map into an MR image for determination of the shape of a tumor in a 
more sensitive imaging modality is also a very promising feature for prostate cancer diagnostic. 

Guidance to predefined targets could increase the sensitivity of prostate biopsies. Suspicious lesions identified in 
MR volumes can be projected into the tracking images and hence enable them to be reached with high precision. If a 
biopsy series needs to be repeated, the clinician can avoid sampling previously examined tissues again, which leads to 
a higher coverage of the gland. A limitation of the system is, however, that it is not yet possible to perform guidance 
to targets in real-time. 

To date, the industrial version of the presented biopsy assistance system has been used on more than 200 biopsy 
sessions and has proven to work reliably in clinical practice. We are currently putting together the protocols for in- 
depth assessment of the clinical potential of the system. We are particularly interested in the impact of the system 
on the sensitivity of US guided prostate biopsies and biopsy repetitions. A preliminary study [20] showed that it 
is difficult for clinicians to reach the zonal targets of the clinical 12-core standard protocol accurately. Giving the 
cUnician a visual feed-back about the sample distribution after the intervention lead already to a significant learning 
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effect. We are also interested in the sensitivity and specificity of MR-targeted biopsies, which is an interesting option 
when the initial biopsy series was negative while the cancer suspicion could not be discarded. 

In the longer term, it would be interesting to consider the usage of the prostate tracking system for therapy, 
knowing that a number of therapeutic interventions on the prostate are performed under endorectal US control (e.g. 
brachytherapy, HIFU, cryotherapy, ...). Focal therapy has the potential to lead to less side-eflFects than radical therapy, 
but it requires the accurate positioning of the therapeutic instrument to be safe, reliable and eff'ective. US based 
prostate tracking can help to achieve this accuracy. 
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